home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / nurb_polyg / NurbRefine.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  3.4 KB  |  136 lines
  1. /*
  2.  * NurbRefine.c - Given a refined knot vector, add control points to a surface.
  3.  *
  4.  * John Peterson
  5.  */
  6.  
  7. #include <stdlib.h>
  8. #include <stdio.h>
  9. #include <math.h>
  10.  
  11. #include "nurbs.h"
  12.  
  13. /*
  14.  * Given the original knot vector ukv, and a new knotvector vkv, compute
  15.  * the "alpha matrix" used to generate the corresponding new control points.
  16.  * This routines allocates the alpha matrix if it isn't allocated already.
  17.  *
  18.  * This is from Bartels, Beatty & Barsky, p. 407
  19.  */
  20. static void
  21. CalcAlpha( double * ukv, double * wkv, long m, long n, long k, double *** alpha )
  22. {
  23.     register long i, j;
  24.     long brkPoint, r, rm1, last, s;
  25.     double omega;
  26.     double aval[MAXORDER];
  27.  
  28.     if (! *alpha)    /* Must allocate alpha */
  29.     {
  30.     CHECK( *alpha = (double **) malloc( (long) ((k+1) * sizeof( double * ))) );
  31.     for (i = 0; i <= k; i++)
  32.         CHECK( (*alpha)[i] = (double *) malloc( (long) ((m + n + 1)
  33.                             * sizeof( double ))) );
  34.     }
  35.  
  36.     for (j = 0; j <= m + n; j++)
  37.     {
  38.     brkPoint = FindBreakPoint( wkv[j], ukv, m, k );
  39.     aval[0] = 1.0;
  40.     for (r = 2; r <= k; r++)
  41.     {
  42.         rm1 = r - 1;
  43.         last = MIN( rm1, brkPoint );
  44.         i = brkPoint - last;
  45.         if (last < rm1)
  46.         aval[last] = aval[last] * (wkv[j + r - 1] - ukv[i])
  47.                 / (ukv[i + r - 1] - ukv[i]);
  48.         else
  49.         aval[last] = 0.0;
  50.  
  51.         for (s = last - 1; s >= 0; s-- )
  52.         {
  53.         i++;
  54.         omega = (wkv[j + r - 1] - ukv[i]) / (ukv[i + r - 1] - ukv[i]);
  55.         aval[s + 1] = aval[s+1] + (1 - omega) * aval[s];
  56.         aval[s] = omega * aval[s];
  57.         }
  58.     }
  59.     last = MIN( k - 1, brkPoint );
  60.     for (i = 0; i <= k; i++)
  61.         (*alpha)[i][j] = 0.0;
  62.     for (s = 0; s <= last; s++)
  63.         (*alpha)[last - s][j] = aval[s];
  64.     }
  65. }
  66.  
  67. /*
  68.  * Apply the alpha matrix computed above to the rows (or columns)
  69.  * of the surface.  If dirflag is true do the U's (row), else do V's (col).
  70.  */
  71. void
  72. RefineSurface( NurbSurface * src, NurbSurface * dest, Boolean dirflag )
  73. {
  74.     register long i, j, out;
  75.     register Point4 * dp, * sp;
  76.     long i1, brkPoint, maxj, maxout;
  77.     register double tmp;
  78.     double ** alpha = NULL;
  79.  
  80.     /* Compute the alpha matrix and indexing variables for the requested direction */
  81.  
  82.     if (dirflag)
  83.     {
  84.     CalcAlpha( src->kvU, dest->kvU, src->numU - 1, dest->numU - src->numU,
  85.            src->orderU, &alpha );
  86.     maxj = dest->numU;
  87.     maxout = src->numV;
  88.     }
  89.     else
  90.     {
  91.     CalcAlpha( src->kvV, dest->kvV, src->numV - 1, dest->numV - src->numV,
  92.            src->orderV, &alpha );
  93.     maxj = dest->numV;
  94.     maxout = dest->numU;
  95.     }
  96.  
  97.     /* Apply the alpha matrix to the original control points, generating new ones */
  98.  
  99.     for (out = 0; out < maxout; out++)
  100.     for (j = 0; j < maxj; j++)
  101.     {
  102.         if (dirflag)
  103.         {
  104.         dp = &(dest->points[out][j]);
  105.         brkPoint = FindBreakPoint( dest->kvU[j], src->kvU,
  106.                        src->numU-1, src->orderU );
  107.         i1 = MAX( brkPoint - src->orderU + 1, 0 );
  108.         sp = &(src->points[out][i1]);
  109.         } else {
  110.         dp = &(dest->points[j][out]);
  111.         brkPoint = FindBreakPoint( dest->kvV[j], src->kvV,
  112.                        src->numV-1, src->orderV );
  113.         i1 = MAX( brkPoint - src->orderV + 1, 0 );
  114.         sp = &(src->points[i1][out]);
  115.         }
  116.         dp->x = 0.0;
  117.         dp->y = 0.0;
  118.         dp->z = 0.0;
  119.         dp->w = 0.0;
  120.         for (i = i1; i <= brkPoint; i++)
  121.         {
  122.         tmp = alpha[i - i1][j];
  123.         sp = (dirflag ? &(src->points[out][i]) : &(src->points[i][out]) );
  124.         dp->x += tmp * sp->x;
  125.         dp->y += tmp * sp->y;
  126.         dp->z += tmp * sp->z;
  127.         dp->w += tmp * sp->w;
  128.         }
  129.     }
  130.  
  131.     /* Free up the alpha matrix */
  132.     for (i = 0; i <= (dirflag ? src->orderU : src->orderV); i++)
  133.     free( alpha[i] );
  134.     free( alpha );
  135. }
  136.